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Direct detection experiments are poised to detect dark matter in the form of weakly in- 
teracting massive particles (WIMPs). The signals expected in these experiments depend 
on the ultra-local WIMP density and velocity distribution. Firstly we review methods 
for modelling the dark matter distribution. We then discuss observational determina- 
tions of the local dark matter density, circular speed and escape speed and the results 
of numerical simulations of Milky Way-like dark matter halos. In each case we highlight 
the uncertainties and assumptions made. We then overview the resulting uncertainties 
in the signals expected in direct detection experiments, specifically the energy, time and 
direction dependence of the event rate. Finally we conclude by discussing techniques for 
handling the astrophysical uncertainties when interpreting data from direct detection 
experiments. 
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1. Introduction 

Weakly Interacting Massive Particles (WIMPs) are a promising dark matter candi- 
date as they are generically produced in the early Universe with roughly the right 
density. Furthermore supersymmetry (SUSY) provides a well-motivated concrete 
WIMP candidate in the form of the lightest neutralino (for reviews see e.g. Refs. [TJ 
[H[3]). WIMPs can be directly detected in the lab via their elastic scattering off target 
nuclei SI. Numerous direct detection experiments are underway and they are prob- 
ing the regions of WIMP mass-cross-section parameter space populated by SUSY 
models (see e.g. Ref. [5]). 

The realisation that uncertainties in the velocity distribution, /(v), will affect 
the direct detection signals dates back to some of the earliest direct detection papers 
written in the 1980s (e.g. Ref. |6]). We first discuss the standard halo model and 
other approaches to modelling the Milky Way halo (Sec. [2}. We then discuss what 
observations and simulations tell us about the dark matter distribution (Sees. |3] 
and 2] respectively) . We then examine the direct detection signals (Sec. [5]) and how 
the astrophysical uncertainties affect these signals (Sec. [6]). We conclude in Sec. [7] 
by discussing strategies for handling the uncertainties, including both 'astrophysics 
independent' approaches and parameterizing the WIMP speed distribution. 
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2. Halo modelling 

2.1. Standard halo model 

The steady-state phase space distribution function, /, of a collection of collisionless 
particles is given by the solution of the collisionless Boltzmann equation: 

In Cartesian coordinates this becomes 

df df <9$d/ „ 

^-+v.^- 7r /=0, 2 
at ax ax av 

where $ is the potential. 

The standard halo model (SHM) is an isotropic, isothermal sphere with density 
profile p(r) oc r~ 2 . In this case the solution of the collisionless Boltzmann equation 
is a so-called Maxwellian velocity distribution, given by 

/(v)=JVexp(-^£), (3) 

where N is a normalisation constant. The isothermal sphere has a flat rotation curve 
at large radii and the velocity dispersion is related to the asymptotic value of the 
circular speed (the speed with which objects on circular orbits orbit the Galactic 
centre) w c ,oo = y/2/3a. It is usually assumed that the rotation curve has already 
reached its asymptotic value at the Solar radius, r — Rq, so that a — v/3/2 v c where 
v c = v c (Ro) is the local circular speed. In the SHM the peak speed vq and the circular 
speed are identical, v c — vq, and these parameters are often used interchangeably. 
However this is not the case in general, for instance for a NFW ^ density profile 
v c = 0.88«o at r = r s /2 (where the scale radius, r«, is the radius at which the 
logarithmic slope of the density profile is equal to -2) ^1 

The density distribution of the SHM is formally infinite and hence the velocity 
distribution extends to infinity too. In reality the Milky Way halo is finite, and 
particles with speeds greater than the escape speed, v csc (r) = ^/2|$(r)|, will not 
be gravitationally bound to the MW. This is often addressed by simply truncating 
the velocity distribution at the measured local escape speed, v csc = v csc (R n ), so 
that /(v) is given by eq. (0) for |v| < u csc and /(v) = for |v| > v csc . This sharp 
truncation is clearly unphysical. An alternative, but still ad hoc, approach is to 
make the truncation smooth: 

/(v) = J JV [ exp (~^) exp (~S 1 )] ' l v l <w - c ' (4) 
I 0, |v| > v csc . 

Another approach, used in the King model/lowered isothermal sphere, is to modify 
the SHM distribution function so that it becomes zero for large energies 

The standard parameter values used for the SHM are a local density po = 
p(Rq) = 0.3 GeVcm~ 3 , a local circular speed v c = 220 km s -1 ^^, and a local escape 
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speed v csc = 544 km s 1 HI We will discuss the determination of these parameters 
and the associated uncertainties in Sec. [3] 

2.2. Non-standard halo models 

For spherical, isotropic systems there is a unique relationship between the density 
profile and the distribution function known as the Eddington equation (see e.g. 
Ref.il [H]): 
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(5) 
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where vE'(r) = — $(r) + $(r = oo) and e = — -Ekin + ^{r), where E^in is the kinetic 
energy. The Eddington equation has been used to calculate th e speed distribution for 
a range of spherically symmetric density profiles [ 12 | 13 | 14 | 1 5] Ref. [T5l found, using a 
Bayesian analysis incorporating various dynamical constraints, speed distributions 
that are not too dissimilar to the standard Maxwellian (when the same value of 
the circular speed is used). Ref. [16] used the Eddington equation to motivate a 
phenomenological form for /(|v|) which matches the high speed tail of the speed 
distributions found in simulated halos (see Sec. [4}. 

The Osipkov-Merritt models assume that the distribution function depends 
on the energy, e, and angular momentum, L, only via a single parameter, Q, 

where the constant r a is the anisotropy radius. In this model the velocity anisotropy 

where crg^^ are the velocity dispersions in the r, 9 and <fi directions, has a particular, 
radially dependent, form: 

and the distribution function can be calculated from a modified version of the 
Eddington equation (see e.g. Refs. [9j IT2|) . In these models the peak of the speed 
distribution becomes narrower and there is an excess of high speed particles (e.g. 
Ref. [18]). 

In general there is no unique relationship between the density profile and the ve- 
locity distribution for triaxial and/or anisotropic systems. In these cases a common 
approach is to use the Jeans equations to calculate the lower order moments of the 
distribution function. Multiplying the collisionless Boltzmann equation, eq. ((2|), by 
one of the velocity components and integrating gives the Jeans equations, which in 
Cartesian coordinates take the form 

d(pvj) + d{pvjm) + d$ = Q 
dt dxi dxj 
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We now have three equations (corresponding to j = 1,2,3) for six unknowns 
(vf, u|, u|, v±V2, V2V3, V1V3). To make further progress it is therefore necessary to 
make assumptions about the alignment of the velocity ellipsoid, i.e. to choose coor- 
dinates such that ViVj = for i 7^ j. The velocity distribution is then approximated 
by a multivariate gaussian in these coordinates: 



where of = vf. 

Evans, Carollo and de Zeeuw studied the logarithmic ellipsoidal model, the 
simplest triaxial generalisation of the isothermal sphere. They argue that aligning 
the velocity ellipsoid with conical coordinates is physically well motivated. In the 
planes of the halo, conical coordinates are locally equivalent to cylindrical polar 
coordinates. Hence, in this model, the velocity distribution can be approximated by 
a multi-variate gaussian in cylindrical polar coordinates, with velocity dispersions 
depending on the shape of the halo, the velocity anisotropy and also location within 
the halo. Both the width of the speed distribution and the peak speed can change 



Analytic models inevitably make assumptions (regarding the shape and 
anisotropy of the halo, and their radial dependence) which are almost certainly 
not completely valid. The relationship between the observed properties of dark 
matter halos and the velocity distribution is non-trivial; models with the same bulk 
properties (e.g. shape and local velocity anisotropy) can have velocity distributions 
with very different forms. Furthermore all analytic approaches to halo modelling 
rely on the assumption that the phase space distribution function has reached a 
steady state. To some extent analytic halo models have been superseded by results 
from numerical simulations. It should be emphasized (see Sec. 14.21 and 14. 1| that 
velocity distributions from numerical simulations also involve approximations and 
extrapolations. 

3. Observations 
3.1. Density 

The event rate is directly proportional to the local density, pq. The standard value 
used is po = 0.3GeVcm~ 3 = 0.008M©pc~ 3 = 5 x 10~ 25 gcm~ 3 . As discussed 
in Ref. [20l the origin of the use of this particular value, rather than say po = 
0.4GeVcm~ 3 , is unclear. 

The local density is calculated via mass modelling of the Milky Way (MW). 
This involves constructing a model of the MW (including its luminous components) 
and then finding the range of values of the local density that are consistent with all 
observational data (including, for instance, rotation curve meas uremen ts , velocity 
dispersions of halo stars, local surface mass density, total mass) I 21 | 22 | As empha- 
sised in Ref. [22] the shape of the MW halo has a crucial effect on the local density 




(10) 



(e.g. Refs.HSHHD. 



Astrophysical uncertainties on direct detection experiments 5 



extracted. For a fixed circular speed, in a flattened halo the local density, in the 
Galactic plane, is higher. 

This work has recently been updated to include new data sets, models for the 
MW halo motivated by numerical simulations, and, in some cases, Bayesian sta- 
tistical techniques. Ref. .23, assuming a spherical halo with a cuspy density profile 
(p(r) oc r~ a as r —> 0) found po = 0.30 ± .05 Ge Vcm~ 3 , while Ref. HJJ assuming 
spherically symmetric NFW - and Einasto | 25|26 | p ro gj es f or the MW halo found, 
po = 0.39 ± 0.03GeVcm~ 3 . While these determinations have ~ 10% errors, they 
differ from each other by significantly more than this. This suggests that the sys- 
tematic errors are bigger than the statistical errors. Indeed Ref. [27| finds, considering 
a range of density profiles including cored profiles (p(r) ~ const as r — > 0), values 
in the range po — 0.2 — 0.4 GeVcm~ 3 . 

Other recent work has investigated 'model independent' techniques, which don't 
involve global mass-modelling of the galaxy, and hence have smaller hidden sys- 
tematic errors. Ref. [20] proposed using the equation of centrifugal equilibrium and 
subtracting the contribution of the stellar component of the MW. The resulting 
determination of the local density is po = 0.43 ± 0.11 ± O.lOGeVcm -3 , where 
the uncertainties come from the uncertainty in the slope of the circular speed 
at the Solar radius and the ratio of the Solar radius to the length scale of the 
thin stellar disk. Ref. [28] introduced a method which involves solving the Jeans 
and Poisson equations with minimal assumptions. Using Hipparcos data they find 
po = 0.1llol27 GeVcm -3 , at 90% confidence, assuming the stellar tracer popula- 
tions are isothermal and po = 1.25^0 34 GeV cm -3 if they have a non-isothermal 
profile. An accurate determination of the vertical dispersion profile of the tracer 
population is therefore required to make an unbiased estimate of the local density 
using this method. 

Ref. 29 used a high-resolution simulation of a Milky Way like galaxy, including 
baryons to investigate the effect of halo shape on the local density. Specifically 
they examined how the local density varies from the density averaged in a spherical 
shell, as determined by observations. They find that the density within the stellar 
disk, at a distance i?o = 8kpc from the centre, is a factor of between 1.01 and 1.41, 
with an average of 1.21, larger than the value averaged in a spherical shell. 



3.2. Circular speed 

The local circular speed, v c = v c (Rq), is an important quantity to determine. It 
appears in the Galilean transformation of the velocity distribution into the lab frame 
(see Sec. [5]). It is also related to the radial component of the velocity dispersion, oy, 
by one of the Jeans equations (e.g. Ref. [S]): 

p dr t r ' 

where the anisotropy parameter j3 is defined in eq. ([7]). The SHM has p(r) cx r~ 2 
and the the velocity distribution is isotropic (oy = og = a^, so that j3 = 0) and 
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independent of radius, so that cr r = v c /\/2. 

The standard value of v c of 220kms~ dates back to a 1980s review of the 
Galactic constants and was found by taking an average of the values found from 
a wide range of different analyses. Note that the ratio v c /Ro is better determined 
than either v c or R individually (e.g. Refs. I3T | l32 |) . 

A recent analysis using measurements of the m otio ns of Galactic masers, found 
a significantly higher value, v c — (254 ± 16) kms -1 ,^]. Ref. [33] reanalyzed the data 
using a more general model for the maser velocity distribution (including allowing 
a non-zero velocity dispersion tensor) and found that the maser data places only a 
relatively weak constraint on v c . When combined with other measurements (from 
the proper motion of Sgr A*, and the orbit of the GD-1 stellar stream), they found 
v c = (236 ± lljkms" 1 , assuming a flat rotation curve. Meanwhile Ref. 1321 found the 
value of v c determined from the maser data depends strongly on the MW model 
used. Using a range of models for the rotation curve, including a power-law with free 
slope, they found values in the range v c = (200±20)kms~ 1 to v c = (279±33)kms~ 1 . 
This illustrates that, as in the case of the local density, systematic, modelling errors 
are important. 



3.3. Escape speed 

The escape speed is the speed required to escape the local gravitational field of 
the MW, v csc (r) = y/2\$(r)\. The local escape speed, v csc = v csc (R ), is estimated 
from the speeds of high velocity stars. To do this it is necessary to parameterise 
the shape of the high speed tail of the velocity distribution. Assuming that the 
velocities are isotropic and the Jeans theorem applies, the asymptotic form of the 
velocity distribution can be written as I 34 | 35 | 

n J(«c 2 S c-|v| 2 ) fc = [(^c-|v|)K sc + |v|)] fc , |v|< Uo8C) 

/(M)cx< (12) 
[0, |v| > v csc . 

Traditionally a value i> OS c = 650kms~ 1 , corresponding to the upper 90% confidence 
limit from Ref. [3H has been used. 

Ref. im has updated these measurements, using additional data from the RAVE 
survey and using a prior on k, k <G [2.7, 4.7] (motivated by analysis of the speed 
distributions of stellar particles in simulated halos) . They find that the escape speed 
lies in the range 498 km s" 1 < t> csc < 608 km s -1 at 90% confidence, with a median 
likelihood of v csc = 544krus~ . 



4. Simulations 

A number of high resolution, dark matter only, simulations of the formation of M ilky 
Way-like halos in a cosmolo gica l context have been carried out (e.g. Aquarius I^SI 
GHALO ^3 and Via Lactea^^j). The veloci ty distribut ions of these halos deviate 
systematically from a multivariate Gaussian I 39 | 40 | 4 1]8| There are more low speed 
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particles and the peak in the distribution is lower and natter in shape (i.e. the 
distribution is platykurtic)H. Several fitting lunctions have been considered; a Tsallis 
distribution (which arises from non-extensive statistical mechanics) ^ and a 
modified Maxwellian HOUI 



/Or) = exp 
/K) = — exp 
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Va + vi are the radial and tangential components of the veloc- 



where v r and Vt = 

ity, N T u are normalisation factors and a Tt t are free-parameters which parameterise 
the deviations from a Maxwellian distribution. The most likely speed deviates from 
the circular speed: for VL2 and GHALO v c /vq s» 0.85 and 0.86 respectively^. 

The velocity distributions also have stochastic features at high speeds. There are 
broad bumps which vary from halo to halo, but are independent of pos ition within a 
given halo and are thought to reflect the formation history of the halo I4US1 Kuhlen 
et al. ^ also find narrow spikes in some locations, corresponding to tidal streams. 

Ref. [16] presented an ansatz for the velocity distribution 
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(15) 



which fits the high speed tail of the smooth component of the speed distributions 
found in simulations. The shape of the high speed tail of the distribution is de- 
termined by the slope of the density profile at large radii. Using the Eddington 
equation the parameter k can be related to the outer sl ope of the density profile, 
7, (p{r) oc r~ 7 for large r), by k = 7 — 3/2 for 7 > 3 As 7 — > 3 (the value 
corresponding to the NFW profile^) the calculation breaks down, and numerical 
fits to eq. (fT~5|) find k ~ 2. More generally for outer slopes in the range 7 ~ 3 — 5, k 
lies in the range k = [1.5 — 3.5]. Eq. ([i"5"j) with k in this range provides a good fit to 
the high speed tails of the speed distributions found in simulation and has less high 
speed particles than the tail of the standard Maxwellian distribution with a sharp 
cut-off. Note, however, that eq. <p~5]) does not match the low and moderate speed 



behaviour of the simulation speed distributions, possibly due to the assumption of 



isotropy 
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4.1. Effects of baryons 

The simulations discussed above contain dark matter only, while baryons dominate 
in the inner regions of the Milky Way. Simulating baryonic physics is extremely dif- 
ficult, and producing galaxies whose detailed properties match those of real galaxies 

a Note, however, that the deviation is smaller in the lab frame than in the Milky Way rest frameEI. 
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is a major challenge. Some recent simulations have found that late merging sub- 
halos are preferentially dragged towards the disc by dynamical friction, w here they 
are destroyed leading to the formation of a co-rotating dark disc (DD) | 44 | 45 | 46 [ 
Ref. H7] modelled the DD velocity distribution as a gaussian with isotropic disper- 
sion, (Tdd = 50kms _1 and lag v\ ag = 50kms _1 , matching (roughly) the kinematics 
of the Milky Way's stellar thick disc and considered DD densities in the range 
0.5 < (pdd/ph) < 2, where pu is the local halo density. 

The properties (and even existence) of the dark disc are highly uncertain. Ref. 1481 
argues that to be consistent with the observed morphological and kinematic proper- 
ties of the Milky Way's thick disc, the Milky Way's merger history must be quiescent 
compared with typical ACDM merger histories, and hence the DD density must be 
relatively small, pdd < 0.2pn- The total (halo + disc) local density can be probed 
by the kinematics of stars (e.g. Ref. |49)) . There is no evidence for a dark matter disc 
(i.e. the data are well fit without a dark disc) and a thick, dense, (pdd/ph) > 0.5, 
dark disc is excluded Refs. l48l and I5T1 also argue that the DD velocity dispersion 
is likely to be substantially larger than that of the stellar thick disc. 

Baryonic physics will also affect the halo speed distribution. For instance gas 
cooling makes halos more spherical (e.g. Ref.l52l). The local speed distribution found 
in Ref. [46] is well fit by a Tsallis distribution, but appears to deviate less from the 
standard Maxwellian than the distributions found in dark matter only simulations. 



4.2. Ultra-local structure 

A further caution is that the scales resolved by simulations are many orders of 
magnitude larger than those relevant for direct detection experiments. The Earth 
moves at - 200 krns" 1 - 0.1 mpc yr 1 , therefore direct detection experiments probe 
the dark matter on sub mpc scales, while the numerical simulations discussed above 
have gravitational softening of order 10 pc. 

Vogelsberger and White have developed a new technique to follow the fine- 
grained phase sp ace distribution in simulations as it stretches and folds under the 
action of gravity ^3. They find that the median den sity of the resulting streams is 
of order 10~ 7 the local halo density. Schneid er et ah j~ h ave studied the evolution 
of the first, and smallest, roughly Earth mass I 55 | 56 | 57 | 58 | m icrohalos to form in the 
Universe. They find that tidal disruption and encounters with stars produce tidal 
streams with average density ~ 10 -4 the local halo density. These results (see also 
Ref. [59]) suggest that the ultra-local dark matter density and velocity distribution 
should not be drastically different to those on the scales resolved by simulations. The 
ultra-loc al velocity distr ibution may, however, contain some features or fine-grained 
structure ™Ea53EI. 
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5. Differential event rate and signals 

The differential event rate for elastic scattering, assuming spin-independent coupling 
with identical couplings to the proton and neutron, is given by (e.g. Refs. IT1 165|): 



where po is the local WIMP density, /(v, i) the WIMP velocity distribution in 
the lab frame, a p the WIMP scattering cross section on the proton, fi px = 
(m p m x ) I (m p + m x ) the WIMP-proton reduced mass, A and F(E) the mass num- 
ber and form factor of the target nuclei respectively and E is the recoil energy. The 
lower limit of the integral, v m i n , is the minimum WIMP speed that can cause a 
recoil of energy E: 



where uia is the atomic mass of the detector nuclei and /ia x the WIMP-nucleon 
reduced mass. The event rate for inelastic scattering of WIMPs also depends on 
the WIMP density and velocity distribution. However, due to the altered kinematics 
the relationship between v m i n and E is changed, and hence the effects of changes in 
the velocity distribution are different I 67 | 8 | 51 | 6 8] 

The WIMP speed distribution must be transformed from the Galactic rest frame 
to the lab frame. This is done by carrying out, a time dependent, Galilean trans- 
formation: v — > v = v + v e (i) lj- The Earth's motion relative to the Galactic 
rest frame, v e (i), is made up of three components: the motion of the Local Stan- 
dard of Rest (LSR), vlsr, the Sun's peculiar motion with respect to the LSR, v Q , 
and the Earth's orbit about the Sun, v° rb ,. The motion of the LSR is defined as 
v lsr = (0,v c ,0), where v c is the local circular speed (see Sec. 13.21 for a discussion 
of recent determinations). The most recent determination of the Sun's motion with 
respect to the LSR, taking into account the effects of the metallicity gradient in the 
disc, finds v = (Z7 , V e , W Q ) = (11.1, 12.2, 7.3) kms" 1 71 in Galactic co-ordinates 
(where U points towards the Galactic center, V is the direction of Galactic rotation 
and W towards the North Galactic Pole). Note that the value of Vq is significantly 
(~ 7kms~ 1 ) larger than previously found The new larger value is also supported 
by the analysis of the motions of masers EH. Accurate expressions for the Earth's 
orbit can be found in Ref. 1731 Simpler expressions, which are acceptable for most 
practical purposes, can be found in Refs. [7H EH 

The differential event rate, eq. (fT6|) . depends on the target nuclei mass, the (a 
priori unknown) WIMP mass and the integral of the velocity distribution. It is 

b Formally gravitational focusing by the Sun should be taken into account | 69 | 70 | nowever the 
resulting modula tion in the differential event rate is small and only detectable with a very large 
number of events ESI. 




(16) 




(17) 
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therefore useful to rewrite eq. (|16l) as 

dR 



dE 



(E,t) = C x ,AP0 5(Wmin,i) , 



where 



and the prefactor 



/(V,t), 3 



d d v, 



c 



I v| >u n 



A 2 F 2 (£) 



(18) 



(19) 



(20) 



contains the WIMP and target dependent terms and is independent of the astro- 
physical WIMP distribution. 

We will now discuss the energy (Sec. l5.1[ ). time fSec. l5.2[) and direction (Sec. l5.3|) 
dependence of the differential event rate. In order to make concrete statements we 
will for now assume the standard Maxwellian velocity distribution, eq. (J3J. In Sec. [5] 
we will discuss how uncertainties in the WIMP distribution affect these signals. 



5.1. Energy dependence 

The shape of the energy spectrum depends on both the WIMP mass and the mass 
of the target nuclei. This can be seen, following Lewin and Smith IJi = !J , by assuming 
a standard Maxwellian velocity distribution, eq. ([3]), and, initially, neglecting the 
Earth's velocity and the Galactic escape speed. This allows the time averaged energy 
spectrum to be written as 

>=G§)-(-il)^ 

where (dR/dE)o is the event rate in the E — > keV limit, and _Er, the characteristic 
energy scale, is given by 

E R = -^L± . (22) 
m A 

When the Earth's velocity and the Galactic escape speed are taken into account 
eq. ([2~Tj) is still a reasonable approximation to the event rate if (dR/dE)o and Er 
are both multiplied by constants of order unity. If m x <C toa, En, oc m 2 , while if 
m x ?7iA, is independent of the WIMP mass. This indicates that the WIMP 
mass can be determined from the ener gy spectru m, provided it is not significantly 
larger than the mass of the target nucleilZSEZCS] Furthermore measuring consistent 
spectra for two different target nuclei could in principle confirm the WIMP origin of 
these spectra (e.g. Ref. [79}. This is sometimes referred to as the 'materials signal'. 
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5.2. Annual modulation 

Due to the Earth's orbit about the Sun, the net velocity of the lab with respect to 
the Galactic rest frame varies annually. The net speed is largest in the Summer and 
hence there are more high speed WIMPs, and less low spee d WI MPs, in the lab 
frame. This produces an annual modulation of the event rate ISlSOl The differential 
event ra te p eaks in Winter for small recoil energies and in Summer for larger recoil 
energies I^JJ The energy at which the annual modulation changes phase is often 
referred to as the 'crossing energy'. Its value depends on the WIMP and target 
masses, and hence could be used to determine the WIMP mass ^21. 

Since the Earth's orbital speed is significantly smaller than the Sun's circular 
speed the amplitude of the modulation is small and, to a first approximation, the 
differential event rate can, for the standard halo model, be written approximately 
as a Taylor series: 

^(E,t)^^(E)[l + A(E)cosa(t)} , (23) 

where a(t) = 2ir(t — to)/T, T = 1 year and to ~ 150 days. For the standard halo 
model the fractional amplitude of the modulation is approximately given by 



A(E) 



-0.034 1 - 



^0-014 -1)^ + 3.7) x p <x<z 



where x = V m i a /v c , x p = 0.89 is the value of x at which the sign of the modu- 
lation reverses and z = v csc /v c . For u m i n > u C sc, in the extreme tail of the speed 
distribution, the shape of the modulation is non-sinusoidal. 

For very small energies A(E) is negative (i.e. the maximum occurs in December 
rather than June). As E is increased, |A(i£)| initially decreases to zero at which 
point the phase of the annual modulation changes so that the maximum occurs in 
Summer. As E is increased further the fractional amplitude continues increasing. 
This is potentially misleading however, as the mean event rate, and hence the raw 
amplitude, becomes very small at large E. After the phase change on increasing E 
further the raw amplitude increases to a local maximum before decreasing again 
and tending to zero. For the SHM, for measurable energies, the amplitude of the 
modulation lies in the range 1 — 10%. 



5.3. Direction dependence 

Our motion with respect to the Galactic rest frame also produces a direction depen- 
dence of the event rate. The WIMP flux in the lab frame is peaked in the direction 
of motion of the Sun (towards the constellation CYGNUS), and hence the recoil 
spectrum is, for most energies El, peaked in the direction opposite to this EH. This 



c At low en ergies the maximum recoil rate is in a ring around the average WIMP arrival direc- 
tionEESl 
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directional signal is far larger than the annual modulation; the event rate in t he 
backward direction is several times larger than that in the forward direction A 
detector which can measure the recoil directions is required to detect this signal 
(see Ref. [53 for an overview of the status of directional detection experiments) . 

The full direction dependence of the event rate is most co mpa ctly written in 
terms of the radon transform of the WIMP velocity distribution 1221 

dR po<r P A 2 



dE dn 4tt/x 



-F 2 (E)f(v min ,q), (24) 



px'"x 

where dil = d0dcos7, q is the recoil direction and /(f m in,q) is the 3-dimensional 
Radon transform of the WIMP velocity distribution 

/(umin,q) = J <S(v.q-v min )/(v)d 3 v. (25) 

Geometrically the Radon transform, /(w m ; n ,q), is the integral of the function /(v) 
on a plane orthogonal to the direction q at a distance v m in from the origin. See 
Ref. 89 for an alternative, but equivalent, expression. 

While the directional recoil rate depends on both of the angles which specify a 
given direction, the strongest signal is the differential of the event ra te with respect 
to the angle between the recoil and the direction of solar motion, 7, 1 86 | 90 | 

d2R P0<Tp A*F*(E) I'* f(v min A)^. (26) 



dEd COS7 47T/i px //(, x J 

For the standard Maxwellian velocity distribution | 86 | Q0 | 



d 2 R C(x,A)p 
- — = 1= exp 

d£/dcos7 V7rt; c 



(27) 



where u° rb 'P is the component of the Earth's velocity parallel to the direction of 
Solar motion. 

An ideal detector capable of measuring the nuclear recoil vectors (including 
their senses, +q versus — q) in 3-dimensions, with good angu lar resolu tion, could 
reject isotropy of WIMP recoils with only of order 10 events I 91 | 89 | 92 | ]yj os ^ ; \y U i 
not all, backgrounds would produce an isotropic Galactic recoil distribution (due 
to the complicated motion of the Earth with respect to the Galactic rest frame). 
An anisotropic Galactic recoil distribution would therefore provide strong, but not 
conclusive, evidence for a Galactic origin of the recoils. Roughly 30 events would 
be required for an ideal detector to confirm that the peak recoil direction coincides 
with the inverse of the directio n of Solar motion, hence confirming the Galactic 
origin of the recoil events I 93 | 94 [ 



6. Consequences 

In this section we will discuss how the uncertainties described in Sees. |3] and [4] affect 
the signals expected in direct detection experiments. 
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6.1. Energy dependence 

Since the normalisation of the event rate is directly proportional to the product of 
the cross-section, a p , and the local density, po, the uncertainty in po leads directly 
to an uncertainty in measurements of, or constraints on, a p . Since the density is 
not expected to vary on mpc scales, the uncertainty in er p is the same for all direct 
detection experiments, and does not affect comparisons of e.g. the exclusion limits 
from different experiments. It does however affect comparisons with collider and 
other constraint s on a p . Furthermore it can significantly bias the determination of 
the WIMP mass^. 

The characteristic energy scale, Er, given by eq. (|2"2"1) . depends on the local 
circular speed. Therefore the uncertainty in the local circular speed, v c , lea ds to a 
0(10%) uncertainty in the differential event rate, and hence exclusion limits EEl2fil 
The nature of the c hange in the exclusion limits is similar, but not identical, for 
differ ent exp eriments ^ ! ^ . It will also lead to a bias in determinations of the WIMP 
mass I 76|77 | This can be seen by differentiating the expression for the characteristic 
energy, eq. 



1 



TOA 



^ . (28) 



Since the differential event rate is given by an integral over the velocity dis- 
tribution the differentia l event rate depends only weakly on the detailed form of 
the v elocity distribution I 97 | 98 [ Consequently the resulting uncertainty in exclusion 
limits Q2168I an( j determinations of the WIMP mass is usually fairly small. There 
are, however, some exceptions to this statement. Firstly if the WIMP is light and/or 
the experimental energy threshold is high compared with the characteristic energy, 
then only the high speed tail of the speed distribution is probed. In this case the 
uncertainties in its shape can have a significant effect on the expected energy spec- 
trum and hence exclusions limits or allowed regions I 68 | 16 [ Secondly if there is a 
dark disc there will be an additional population of low speed WIMPs. If the dark 
disc density is sufficiently high and its velocity dispersion sufficiently small this will 
significantly change the energy spectrum and hence exclusion limits and mass de- 
terminations The size of these change s depends on the properties of the dark 
disc, which are currently uncertain I 48 | 51 | j± w±MP stream, with sufficiently high 
density, would add a sloping step to the differential event rate'-'- 2 '. 



6.2. Annual modulation 

The annual modulation arises from the, relatively small, shift in the speed distribu- 
tion in the lab frame over the course of the year. It is therefore far more sensitive to 
the speed distribution than the time averaged energy spectrum. Both the amplitude 
and phase of the modulation, and hence the regions of parameter space compatible 
with an observed signal, can vary significantly (e.g. Refs. I100[ I10H I102[ 1103} 
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The uncertainty in the value of v c can change the amplitude by a factor of order 
unity, while the uncertainty i n the sh ape of the halo velocity distribution changes 
the amplitude by ~ 0(10%) | 10^ I 96 | A high density dark disc could change the 
annual modulation signal significantly, producing a 2nd, low energy, m axim um in 
the amplitude of the modulation and changing the phase significantly 1221. For a 
WIMP stream the position and height of the step in the energy spectrum produced 
would vary annually (e.g. Ref . |83|) . 



6.3. Direction dependence 

The detailed direction dependence of the e vent rate is sensitive to the velocity dis- 
tribution, however the anisotropy is robust [HUEHIHH The number of events required 
to detect anisotropy or to demonstrate that the median inverse rec oil direction 
coincides with the direction of solar motion only vary by of order 10% | 91 | 89 | 92 |94| 
The features in the high speed tail can cause the median inverse direction of high 
energy recoils to deviate from the direction of solar motion A stream of WIMPs 
would produce recoils which are strongly peaked in the opposite direction l^QI. In 
the long term, studying the median direction of high energy recoils could allow high 
speed features to be detected and the formation history of the Milky Way probed. 



7. How to handle the uncertainties? 
7.1. Astrophysics independent methods 

It is possible to carry out an 'astrophysics independent' comparison of data from 
multiple experiments. This effectively involves making comparisons of the energy 
spectra, so that the integral over the velocity distribution, g(v m in) defined in eq. (fl9|) . 
cancels. For instance, with data from two different experiments using two different 
targets, the WIMP mass can in principle be determ ined , without any assumptions 
about /(v), by taking moments of the energy spectra^l However this method leads 
to a systematic underestimate of the WIMP mass if it is comparable to or larger 
than the mass of the target nuclei. 

Ref. 11081 see also Ref. 11091 showed how the differential event rates in two ex- 
periments, '1' and '2' with target nuclei with mass numbers A\ and are related. 
If the first experiment is sensitive to recoil energies in the range [E l ° w ,E^ E ] it 
probes g(v min ) for v l ^ n < v min < v^f n , where v min is related to the recoil energy, 
WIMP mass and target nuclei mass by eq. (fl7]) . In experiment 2 these values of 
Vmin correspond to recoils in the range 

2 

r^low E U S h , = ^A2 X ™A1 [£low ^highj _ ^ 
H A l x m A2 



For recoil energies in this range the differential event rates in the two experiments 
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are then related by 

¥k(E ) = F " {E2) dRl 

dE ( 2) C x , Al F 2(rt^ E ) dE 

This relation depends on the (unknown) WIMP mass, therefore the comparison has 
to be made separately for a range of m x v alues . This method has recently been 
extended to include annual modulation data HA^. ^^^^^^^ 

This approach has proved particularly useful m comp aring the 

event ra te ex cesses and ann ual modulations found in the CoGeNT | H3 | 114 | 
CRESST ED I and D AMA El experiments with exclusions limits from CDMsQlZl 
and Xenon l 118 l 119 l. Ref. QT0] found that the DAMA, CoGeNT and CRESST-II 
data are compatible which each other, but not CDMS and XENON, if the velocity 
distribution is very anisotropic (which leads to a large modulation fraction). 



( ^L E2 ) . (30) 

\^A2 x m Al 



7.2. Parameterizing the speed distribution 

The astrophysics independent methods are invaluable for assessing the compatibility 
of data from different experiments. However the WIMP mass is not a priori known, 
and the goal of WIMP direct detection experiments is not just to detect WIMPs, 
but to also measure their mass and cross-section. While there are uncertainties in 
the velocity distribution it is not completely unknown, and therefore does not need 
to be completely removed from the analysis. 

Strigari and Trotta^lhave shown how the WIMP properties can be constrained 
by a single future direct detection experiment when astronomical data (such as the 
kinematics of halo stars) are used to jointly constrain the WIMP properties and 
the parameters of a model for the Milky Way halo. They assumed an isotropic 
Maxwellian speed distribution characterised by the peak speed, vq, and the escape 
speed, ? w.. 

Peter | 120 | 121 | nag ex amined how data sets from multiple direct detection ex- 
periments could be used to jointly constrain a Maxwellian parametrisation of the 
WIM P sp eed distribution and the WIMP parameters (mass and cross-section). Pato 
et al. [1221 have taken a similar approach using the parameterization in eq. (|15l) . Pe- 
ter has recently extended this work by allowing the peak speed and circular speed 
to differ (as expected if the density profile is not p(r) tx r -2 ) and also considering 
an empirical speed distribution consisting of a five or ten bin step function 1 I 

Combining data from multiple experiments with different targets significantly 
increases the accuracy with which the WIMP parameters can be measured I 120 | 29 [ 
However fixing the form of the speed distribution leads to biases in the WIMP prop- 
erties, if the true speed distribution differs significantly from the parameterization 
assumed 1 I With a suitable parameterization o the WIMP parameters and speed 
distribution can be jointly probed The form of the optimal parameterization 
of the speed distribution is an open question. 
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8. Summary 

Direct detection event rate calculations often assume the standard halo model, 
with an isotropic Maxwellian velocity distribution with dispersion a = y / 3/2v c = 
270 km s" 1 and a local WIMP density p = 0.3 GeV cm . However it has long been 
realised^ that uncertainties in the speed distribution will affect the signals expected 
in experiments. 

We have discussed the standard halo model and other approaches to halo mod- 
elling, and the assumptions behind them. We then reviewed observational determi- 
nations of quantities that are relevant to direct detection experiments, namely the 
local dark matter density, the local circular speed (which is related to the velocity 
dispersion by the Jeans equations) and the local escape speed. While the statistical 
errors on these quantities are often small the systematic errors, from uncertainties 
in the modelling of the Milky Way, can be significantly larger. 

Next we turned our attention to high resolution, dark matter only simulations 
of the formation of Milky Way-like halos in a cosmological context. They find ve- 
locity distributions that deviate significantly from the standard Maxwellian, and 
have features at high speeds. The effect of baryonic physics on the dark matter 
distribution is not yet well understood. Some recent simulations have found that 
a dark disc may be formed, however its properties (and even existence) are highly 
uncertain. Simulations resolve the dark matter distribution on scales many orders 
of magnitude larger than those probed by direct detection experiments. The latest 
results suggest that the ultra-local dark matter distribution is largely smooth, but 
some features may exist. 

We then reviewed the resulting uncertainties in the energy, time and direction 
dependence of the energy spectrum, and hence constraints on, or measurements of, 
the WIMP parameters. The uncertainty in the local circular speed has a signifi- 
cant effect on the event rate and hence exclusion limits and determinations of the 
WIMP mass. The uncertainty in the local density leads directly to an uncertainty 
in measurements of, or constraints on, the cross-section. Only the high energy tail 
of the energy spectrum is particularly sensitive to the exact shape of the speed dis- 
tribution. The annual modulation is far more sensitive to the shape of the velocity 
distribution; its amplitude and phase can change significantly. The anisotropy of 
the recoil rate is robust to changes in the velocity distribution, however high speed 
features can change the peak direction of high energy recoils, and hence provide a 
way of probing the formation of the Milky Way. 

Finally we discussed techniques for handling the astrophysical uncertainties 
when analysing direct detection data. Astrophysics independent comparisons be- 
tween different experiments can be made using the integral of the velocity distri- 
bution, g(v m i n ). This approach is extremely useful for assessing the compatibility 
of various experiments, but requires the WIMP mass as input. Parameterising the 
WIMP speed distribution, and using data from multiple experiments to jointly 
constrain the WIMP mass and cross-section and speed distribution is a promising 
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approach, but the optimal form for the parameterisation is not yet known. 
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